# Load behavioral data
load('/Users/noah/Desktop/SGE_data_objects/sge_phase12.Rdata')
data$Fq9_SGE2 <- data$Fq9_SGE2[c(1:5012,5144:5791),]
data$Pp10_SGE2 <- data$Pp10_SGE2[c(1:662,797:7334),]
## load functions & intros object
load('/Users/noah/Desktop/SGE_data_objects/behav_functions.RData')
load('/Users/noah/Desktop/SGE_data_objects/datafornoah.RData')
## number of focals per group
for (x in names(data)) intros[which(intros[,"group"]==x),"num.obs"]=length(levels(as.factor((data[[x]][,"obs.no"]))))
# pull out social behaviors and dominance interactions
social<-lapply(data,convert.soc)
elo<-lapply(data,convert.elo)
groom<-lapply(social,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
dominance<-lapply(data,convert.dom)
prox<-lapply(social,function(x){x[which(x[,5]%in%c("PX","LB","ZZ")),]})
intros$finalElo <- NA
for (i in 1:90){
id=tolower(intros[i,"ID"])
group=intros[i,"group"]
intros[i,"finalElo"]=extract_elo(elo[[group]],ID=id)
}
intros<-intros[order(intros$phase,intros$group),]
## Calculate SGE I hierarchy metrics and generate a dataframe containing all metrics
sge.I.hierarchy.metrics <- matrix(nrow = 18, ncol = 6)
count <- 1
for (i in unique(names(dominance))){
tmp <- as.data.frame(dominance[i])
tmp <- tmp[order(tmp[,1]),]
tmp <- subset(tmp, tmp[,3] != tmp[,4])
seqcheck(tmp[,3],tmp[,4],tmp[,1])
elo<-elo.seq(tmp[,3],tmp[,4],tmp[,1])
elo.mat <- creatematrix(elo)
a <- extract_elo(elo)
b <- steeptest(elo.mat, rep = 1000)$Stp
c <- devries(elo.mat)$`h-modified`
d <- dci(elo.mat)
e <- phi(elo.mat)
f <- ttri(elo.mat)$ttri
g <- stab_elo(elo)
sge.I.hierarchy.metrics[count,] <- c(b,c,d,e,f,g)
count <- count + 1
}
sge.I.hierarchy.metrics <- as.data.frame(sge.I.hierarchy.metrics)
sge.I.hierarchy.metrics$phase <- rep(c(1,2), times = c(9,9))
rownames(sge.I.hierarchy.metrics) <- names(dominance)
colnames(sge.I.hierarchy.metrics) <- c('steepness','linearity','DC','phi','tt','stability','phase')
sge.I.hierarchy.metrics$group <- rownames(sge.I.hierarchy.metrics)
Here I’m calculating a variety of metrics that describe the nature of the dominance hierarchies in SGE I. These include steepness, linearity, direction concordance, phi, triangle transitivity, and stability (descriptions of these hierarchy metrics can be found here).
Correlations between hierarchy metrics.
Panel 1 (both phase); Panel 2 (phase 1); panel 3 (phase 2). Triangle transitivity is not included in the plots as all groups had a value of 1. Similarly, all groups but one had a value of 1 for linearity.
# plot correlations between elos, DS, and ordinal ranks for SGE I
ggpairs(as.data.frame(sge.I.hierarchy.metrics[,c(1,3:4,6)]), title = 'Phases 1 + 2')
ggpairs(as.data.frame(sge.I.hierarchy.metrics[1:9,c(1,3:4,6)]), title = 'Phase 1')
ggpairs(as.data.frame(sge.I.hierarchy.metrics[10:18,c(1,3:4,6)]), title = 'Phase 2')
Distribution of hierarchy metrics across groups. Phi is plotted separately as it’s max value is 0.5 rather than 1 as is the case with the other metrics.
metric_melt <- sge.I.hierarchy.metrics[,c(1:3,5:6)]
metric_melt <- reshape2::melt(metric_melt)
phi_melt <- sge.I.hierarchy.metrics[,4]
phi_melt <- reshape2::melt(phi_melt)
p1 <- ggplot(metric_melt, aes(x=as.factor(variable), y = value)) +
geom_boxplot() +
xlab('hierarchy metric')
p2 <- ggplot(phi_melt, aes(y = value)) +
geom_boxplot() +
xlab('phi')
plot_grid(p1,p2)
Here I’m going to show the correlations between all the behavioral variables being used in the study.
# read in behavioral data
sge.I.behavioral.v2 <- read.delim("~/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt")
sge.I.behavioral.v2$rowOrder <- 1:nrow(sge.I.behavioral.v2)
sge.I.behavioral.v2$aggDiff <- sge.I.behavioral.v2$raw.aggGiven.rate - sge.I.behavioral.v2$raw.aggRec.rate
# PC1.agg and PC1.groom
agg_vector <- c('raw.aggRec.contact.rate', 'raw.aggRec.nocontact.rate', 'raw.aggGiven.contact.rate', 'raw.aggGiven.nocontact.rate', 'aggDiff')
tmp1 <- subset(sge.I.behavioral.v2, phase == '1')
tmp2 <- subset(sge.I.behavioral.v2, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,agg_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,agg_vector])), scale = T)
sge.I.behavioral.v2$PC1.agg <- NA
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='1'),'PC1.agg'] <- pca1$x[,1]
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='2'),'PC1.agg'] <- pca2$x[,1]
groom_vector <- c('raw.groomRec.rate', 'raw.groomGiven.rate')
tmp1 <- subset(sge.I.behavioral.v2, phase == '1')
tmp2 <- subset(sge.I.behavioral.v2, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,groom_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,groom_vector])), scale = T)
sge.I.behavioral.v2$PC1.groom <- NA
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='1'),'PC1.groom'] <- pca1$x[,1]
sge.I.behavioral.v2[which(sge.I.behavioral.v2$phase=='2'),'PC1.groom'] <- pca2$x[,1]
# group-center and scale behavioral variables
sge.I.behavioral.v2$overall.agg.received <- sge.I.behavioral.v2$raw.aggRec.contact.rate + sge.I.behavioral.v2$raw.aggRec.nocontact.rate
sge.I.behavioral.v2$overall.agg.given <- sge.I.behavioral.v2$raw.aggGiven.contact.rate + sge.I.behavioral.v2$raw.aggGiven.nocontact.rate
sge.I.behavioral.v2$overall.groom <- sge.I.behavioral.v2$raw.groomRec.rate + sge.I.behavioral.v2$raw.groomGiven.rate
sge.I.behavioral.v2 <- sge.I.behavioral.v2[order(sge.I.behavioral.v2$group),]
group.center <- function(behaveVar){
test <- vector()
for (i in unique(sge.I.behavioral.v2$group)){
tmp <- subset(sge.I.behavioral.v2, group == i)
tmp2 <- scale(tmp[,behaveVar], scale=F, center=T)
test <- rbind(test, tmp2)
}
sge.I.behavioral.v2[,paste('gc.',behaveVar,sep='')] <- as.numeric(scale(test))
return(sge.I.behavioral.v2)
}
behaviorList <- c('overall.groom',
'overall.agg.given',
'overall.agg.received',
'raw.aggRec.contact.rate',
'raw.aggRec.nocontact.rate',
'raw.aggGiven.contact.rate',
'raw.aggGiven.nocontact.rate',
'raw.groomRec.rate',
'raw.groomGiven.rate',
'aggDiff')
for (i in behaviorList){
sge.I.behavioral.v2 <- group.center(i)
}
sge.I.behavioral.v2 <- sge.I.behavioral.v2[order(sge.I.behavioral.v2$rowOrder),]
Elo vs. agg variables
p1 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.aggDiff, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('aggDiff') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
p2 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=PC1.agg, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('PC1.agg') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
p3 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.overall.agg.given, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('Overall agg given') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
p4 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.overall.agg.received, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('Overall agg received') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
plot_grid(p1,p2,p3,p4, ncol = 2)
Elo vs. groom variables
p1 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.raw.groomGiven.rate, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('Overall groom given') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
p2 <- ggplot(sge.I.behavioral.v2,aes(x=elo, y=gc.raw.groomRec.rate, group = phase)) +
geom_point(aes(color=as.factor(phase))) +
xlab('rank (elo)') +
ylab('Overall groom received') +
stat_smooth(method = 'lm', aes(color=as.factor(phase)), se=F, fullrange = T) +
ggtitle('') +
scale_color_manual(values = c("#fb8474","#748ccc"), name = 'Phase') +
theme_minimal()
plot_grid(p1,p2, ncol = 2)
All pairwise comparisons
ggpairs(sge.I.behavioral.v2[,c(5,26:27,31:40)],
upper = list(continuous = wrap("cor", size = 2)),
lower = list(continuous = wrap("points", size = 0.5))) +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
Now we can look at the PVE and loadings for PC1.agg and PC1.groom.
First we’ll look at PC1.agg
## Generate PCA info for PC1.agg variables for SGE I and II both phases
sgeIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt', header = T)
sgeIbehavior$aggDiff <- sgeIbehavior$raw.aggGiven.rate - sgeIbehavior$raw.aggRec.rate
sgeIIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.II.behavioral.v2.txt', header = T)
sgeIIbehavior$aggDiff <- sgeIIbehavior$raw.aggGiven.rate - sgeIIbehavior$raw.aggRec.rate
agg_vector <- c('raw.aggRec.contact.rate', 'raw.aggRec.nocontact.rate', 'raw.aggGiven.contact.rate', 'raw.aggGiven.nocontact.rate', 'aggDiff')
tmp1 <- subset(sgeIbehavior, phase == '1')
tmp2 <- subset(sgeIbehavior, phase == '2')
tmp3 <- subset(sgeIIbehavior, phase == '1')
tmp4 <- subset(sgeIIbehavior, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,agg_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,agg_vector])), scale = T)
pca3 <- prcomp(as.data.frame(scale(tmp3[,agg_vector])), scale = T)
pca4 <- prcomp(as.data.frame(scale(tmp4[,agg_vector])), scale = T)
pov_table <- matrix(nrow = 20, ncol = 2)
pov_table[,1] <- rep(c('I_I','I_II','II_I','II_II'), times=c(5,5,5,5))
pov_table[1:5,2] <- summary(pca1)$importance[2,]
pov_table[6:10,2] <- summary(pca2)$importance[2,]
pov_table[11:15,2] <- summary(pca3)$importance[2,]
pov_table[16:20,2] <- summary(pca4)$importance[2,]
pov_table <- as.data.frame(pov_table)
pov_table$PC <- rep(1:5, times=4)
pov_table <- pov_table[1:10,]
pov_table$V1 <- as.factor(rep(c(1,2), times=c(5,5)))
ggplot(pov_table, aes(x=PC, y = as.numeric(V2), color = V1)) +
geom_line(size=2) +
geom_point(size=2) +
xlab('PC') +
ylab('Proportion of variance explained by PC') +
labs(color = "Phase") +
theme_minimal() +
scale_color_viridis(discrete = T, option = 'E')
I_I <- cbind(reshape2::melt(pca1$rotation), 'I_I')
I_II <- cbind(reshape2::melt(pca2$rotation), 'I_II')
II_I <- cbind(reshape2::melt(pca3$rotation), 'II_I')
II_II <- cbind(reshape2::melt(pca4$rotation), 'II_II')
loading_list <- list(I_I,I_II,II_I,II_II)
loading_table <- data.table::rbindlist(loading_list, use.names = F)
colnames(loading_table) <- c('behavior', 'PC', 'loading', 'exp_phase')
loading_table <- loading_table[1:50,]
loading_table$exp_phase <- as.factor(rep(c(1,2), times = c(25,25)))
ggplot(loading_table, aes(PC, abs(loading))) +
geom_linerange(aes(x = PC, ymin = 0, ymax = abs(loading), group = behavior), color = "lightgray", size = 1.5, position = position_dodge(0.3)) +
geom_point(aes(color = behavior),position = position_dodge(0.3), size = 3) +
facet_grid(rows = vars(exp_phase)) +
theme_minimal() +
scale_color_viridis(discrete = T, option = 'E')
Now we can look at PC1.groom
## Generate PCA info for PC1.agg variables for SGE I and II both phases
sgeIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.I.behavioral.v2.txt', header = T)
sgeIbehavior$aggDiff <- sgeIbehavior$raw.aggGiven.rate - sgeIbehavior$raw.aggRec.rate
sgeIIbehavior <- read.table('/Users/noah/Desktop/SGE_data_objects/sge.II.behavioral.v2.txt', header = T)
sgeIIbehavior$aggDiff <- sgeIIbehavior$raw.aggGiven.rate - sgeIIbehavior$raw.aggRec.rate
groom_vector <- c('raw.groomRec.rate', 'raw.groomGiven.rate')
tmp1 <- subset(sgeIbehavior, phase == '1')
tmp2 <- subset(sgeIbehavior, phase == '2')
tmp3 <- subset(sgeIIbehavior, phase == '1')
tmp4 <- subset(sgeIIbehavior, phase == '2')
pca1 <- prcomp(as.data.frame(scale(tmp1[,groom_vector])), scale = T)
pca2 <- prcomp(as.data.frame(scale(tmp2[,groom_vector])), scale = T)
pca3 <- prcomp(as.data.frame(scale(tmp3[,groom_vector])), scale = T)
pca4 <- prcomp(as.data.frame(scale(tmp4[,groom_vector])), scale = T)
pov_table <- matrix(nrow = 8, ncol = 2)
pov_table[,1] <- rep(c('I_I','I_II','II_I','II_II'), times=c(2,2,2,2))
pov_table[1:2,2] <- summary(pca1)$importance[2,]
pov_table[3:4,2] <- summary(pca2)$importance[2,]
pov_table[5:6,2] <- summary(pca3)$importance[2,]
pov_table[7:8,2] <- summary(pca4)$importance[2,]
pov_table <- as.data.frame(pov_table)
pov_table$PC <- rep(1:2, times=4)
pov_table <- pov_table[1:4,]
pov_table$V1 <- as.factor(rep(c(1,2), times = c(2,2)))
ggplot(pov_table, aes(x=PC, y = as.numeric(V2), color = V1)) +
geom_line(size=2) +
geom_point(size=2) +
xlab('PC') +
ylab('Proportion of variance explained by PC') +
labs(color = "Phase") +
theme_minimal() +
scale_color_viridis(discrete = T, option = 'E')
I_I <- cbind(reshape2::melt(pca1$rotation), 'I_I')
I_II <- cbind(reshape2::melt(pca2$rotation), 'I_II')
II_I <- cbind(reshape2::melt(pca3$rotation), 'II_I')
II_II <- cbind(reshape2::melt(pca4$rotation), 'II_II')
loading_list <- list(I_I,I_II,II_I,II_II)
loading_table <- data.table::rbindlist(loading_list, use.names = F)
colnames(loading_table) <- c('behavior', 'PC', 'loading', 'exp_phase')
loading_table <- loading_table[1:8,]
loading_table$exp_phase <- as.factor(rep(c(1,2), times = c(4,4)))
ggplot(loading_table, aes(PC, abs(loading))) +
geom_linerange(aes(x = PC, ymin = 0, ymax = abs(loading), group = behavior), color = "lightgray", size = 1.5, position = position_dodge(0.3)) +
geom_point(aes(color = behavior),position = position_dodge(0.3), size = 3) +
facet_grid(rows = vars(exp_phase)) +
theme_minimal() +
scale_color_viridis(discrete = T, option = 'E')
As you can see the loadings are not meaningful with only two variables in the PCA, but the first PC explains ~60% of the variance. Not sure is PC1.groom is a valid composite variable.
One metric we will use to assess the data quality is to split the data and see how correlated the behavioral variables are between the two subsets, where we’d expect the two subsets to be well correlated if we’re capturing ‘true’ behavioral rates.
Because the subsets need to maintain the row structure in order to calculate grooming times, the data were split by obs, alternating across the study duration so both subsets contain data from across the entire study.
##caclulate grooming bouts
## Here I'm going to 1) split the data in half, 2) regenerate groom objects, 3) calculate grooming rates for each half, and 4) see how correlated they are
## Can't split by taking alternating rows becasue the row stricture needs to be maintained in order to calculate mins of grooming, so need to split by something that's unrelated to grooming, I'll sort alphabetically, split, and the resort to original order in each subset (by row order)
data_row <- lapply(names(data),function(x){data[[x]]$rowOrder <- 1:nrow(data[[x]]); return(data[[x]])})
names(data_row) <- names(data)
data_half1 <- lapply(names(data_row),function(x){evens <- seq(2,max(unique(data[[x]]$obs.no)),2);return(data_row[[x]][which(data_row[[x]][,3] %in% evens),])})
names(data_half1) <- names(data_row)
data_half2 <- lapply(names(data_row),function(x){odds <- seq(1,max(unique(data[[x]]$obs.no)),2);return(data_row[[x]][which(data_row[[x]][,3] %in% odds),])})
names(data_half2) <- names(data_row)
social_1 <- lapply(data_half1,convert.soc)
groom_1 <- lapply(social_1,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
names(groom_1) <- names(data)
social_2 <- lapply(data_half2,convert.soc)
groom_2 <- lapply(social_2,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
names(groom_2) <- names(data)
groom.function <- function(groomDF){
groom.tot1<-data.frame(matrix(ncol=19,nrow=360))
colnames(groom.tot1)<-c("ID1","ID2","group","no.obs","bouts","ID1.groom","ID2.groom","min","ID1.groom.min","ID2.groom.min","ID1.rank","ID2.rank","proximity","contact.min","CSI","elo.ID1","elo.ID2","ID1.agg","ID2.agg")
groom.tot1[,"ID1"]<-rep(intros[,"ID"],4)
for (i in 1:180) {groom.tot1[i,"group"]<-intros[which(intros[1:45,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
for (i in 181:360) {groom.tot1[i,"group"]<-intros[46:90,][which(intros[46:90,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
groom.tot1<-groom.tot1[order(groom.tot1[,3],groom.tot1[,1]),];rownames(groom.tot1)<-NULL
for (i in unique(intros[,4])) {intros[1:45,][which(intros[1:45,"ID"]==i),"group"]->grp;
groom.tot1[grep("SGE1",groom.tot1$group),][which(groom.tot1[grep("SGE1",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros[1:45,][which(intros[1:45,"group"]==grp),"ID"],i)}
for (i in unique(intros[,4])) {intros[46:90,][which(intros[46:90,"ID"]==i),"group"]->grp;
groom.tot1[grep("SGE2",groom.tot1$group),][which(groom.tot1[grep("SGE2",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros[46:90,][which(intros[46:90,"group"]==grp),"ID"],i)}
groom.tot<-groom.tot1; rm(groom.tot1)
groom.tot[,1]<-as.character(groom.tot[,1]); groom.tot[,2]<-as.character(groom.tot[,2]); groom.tot[,3]<-as.character(groom.tot[,3]); rm(grp); rm(i)
##caclulate grooming bouts
groom.tot$ID1.intro.no<-NA; groom.tot$ID2.intro.no<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
b1=sum(g[,5]=="GM"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="GM"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,5]=b1+b2;groom.tot[id,6]=b1;groom.tot[id,7]=b2
groom.tot[id,4]=intros[which(intros[,"ID"]==id2&intros$group==grp),"num.obs"]
groom.tot[id,"ID1.rank"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"finalElo"]
groom.tot[id,"ID2.rank"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"finalElo"]
groom.tot[id,"ID1.intro.no"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"intro.no"]
groom.tot[id,"ID2.intro.no"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"intro.no"]
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
groom.tot$ID1.intro.date<-NA; groom.tot$ID2.intro.date<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
groom.tot[id,"ID1.intro.date"]=intros[which(intros[,"ID"]==id1&intros$group==grp),"intro.date"]
groom.tot[id,"ID2.intro.date"]=intros[which(intros[,"ID"]==id2&intros$group==grp),"intro.date"]
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
## calculate who ended the grooming bout
groom.tot$ID1.end.bout<-NA; groom.tot$ID2.end.bout<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
b1=sum(g[,5]=="G-"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="G-"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,"ID1.end.bout"]=b1
groom.tot[id,"ID2.end.bout"]=b2
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
groom.tot$ID1.approach<-NA
groom.tot$ID2.approach<-NA
##caclulate approaches
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-prox[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
b1=sum(g[,5]=="PX"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="PX"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,"ID1.approach"]=b1;groom.tot[id,"ID2.approach"]=b2
}; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
## calculate grooming times
groom.tot$ID1.groom.min<-NA
groom.tot$ID2.groom.min<-NA
groom.tot$min<-NA
for (i in 1:nrow(groom.tot)){
id1<-groom.tot[i,"ID1"]; id2<-groom.tot[i,"ID2"]; t1=0;t2=0
g<-groomDF[[as.character(groom.tot[i,"group"])]]
if (id1%in%c(g[,3],g[,4])&id2%in%c(g[,3],g[,4])){
g<-g[which((g[,4]==id1 & g[,3]==id2 )|(g[,3]==id1 & g[,4]==id2)|g[,3]=="ZZ"),]
for (r in 1:nrow(g)){
if (g[r,5]=="GM" & g[r,3]==id1 & g[r,4]==id2) t1=t1+(g[r+1,2]-g[r,2])
if (g[r,5]=="GM" & g[r,4]==id1 & g[r,3]==id2) t2=t2+(g[r+1,2]-g[r,2])
}
groom.tot[i,"ID1.groom.min"]=t1
groom.tot[i,"ID2.groom.min"]=t2
groom.tot[i,"min"]=t1+t2}
}
rm(t1);rm(t2);rm(i);rm(id1);rm(id2);rm(r);rm(g)
###########
phase <- groom.tot$group
phase <- substr(phase,(nchar(phase)+1)-4,nchar(phase))
groom.tot$phase <- phase
## generate a data frame with minutes of grooming (total), elo, # of obs, and group, per individual
## phase I
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
tmp <- subset(groom.tot, groom.tot$phase == 'SGE1')
tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
tmp.minutes <- cbind(tmp.minutes,tmp2)
tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
tmp.names <- cbind(tmp.names,i)
tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}
sge.I.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.I.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.I.groom.rate.df$elo <- as.numeric(as.character(sge.I.groom.rate.df$elo))
sge.I.groom.rate.df$min <- as.numeric(as.character(sge.I.groom.rate.df$min))
sge.I.groom.rate.df$min.given <- as.numeric(as.character(sge.I.groom.rate.df$min.given))
sge.I.groom.rate.df$min.received <- as.numeric(as.character(sge.I.groom.rate.df$min.received))
sge.I.groom.rate.df$no.obs <- as.numeric(as.character(sge.I.groom.rate.df$no.obs))
sge.I.groom.rate.df$phase <- 1
## now for phase II
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
tmp <- subset(groom.tot, groom.tot$phase == 'SGE2')
tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
tmp.minutes <- cbind(tmp.minutes,tmp2)
tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
tmp.names <- cbind(tmp.names,i)
tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}
sge.II.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.II.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.II.groom.rate.df$elo <- as.numeric(as.character(sge.II.groom.rate.df$elo))
sge.II.groom.rate.df$min <- as.numeric(as.character(sge.II.groom.rate.df$min))
sge.II.groom.rate.df$min.given <- as.numeric(as.character(sge.II.groom.rate.df$min.given))
sge.II.groom.rate.df$min.received <- as.numeric(as.character(sge.II.groom.rate.df$min.received))
sge.II.groom.rate.df$no.obs <- as.numeric(as.character(sge.II.groom.rate.df$no.obs))
sge.II.groom.rate.df$phase <- 2
sge.groom.rate.df <- rbind(sge.I.groom.rate.df,sge.II.groom.rate.df)
return(sge.groom.rate.df)
}
groom.df.1 <- groom.function(groom_1)
groom.df.2 <- groom.function(groom_2)
# turn minutes into rates
groom.df.1$min.rate <- (groom.df.1$min / groom.df.1$no.obs) * 2
groom.df.1$min.given.rate <- (groom.df.1$min.given / groom.df.1$no.obs) * 2
groom.df.1$min.received.rate <- (groom.df.1$min.received / groom.df.1$no.obs) * 2
groom.df.2$min.rate <- (groom.df.2$min / groom.df.2$no.obs) * 2
groom.df.2$min.given.rate <- (groom.df.2$min.given / groom.df.2$no.obs) * 2
groom.df.2$min.received.rate <- (groom.df.2$min.received / groom.df.2$no.obs) * 2
tmp <- cor.test(groom.df.1$min.rate,groom.df.2$min.rate)
p1 <- ggplot(groom.df.1, aes(x=min.rate,y=groom.df.2$min.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
xlab('grooming rate subset 1') +
ylab('grooming rate subset 2') +
ggtitle('total grooming rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
#stat_poly_eq(formula = y~x,
# aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
# parse = TRUE)
tmp <- cor.test(groom.df.1$min.given.rate,groom.df.2$min.given.rate)
p2 <- ggplot(groom.df.1, aes(x=min.given.rate,y=groom.df.2$min.given.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
xlab('grooming rate subset 1') +
ylab('grooming rate subset 2') +
ggtitle('grooming given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
#stat_poly_eq(formula = y~x,
# aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
# parse = TRUE)
tmp <- cor.test(groom.df.1$min.received.rate,groom.df.2$min.received.rate)
p3 <- ggplot(groom.df.1, aes(x=min.received.rate,y=groom.df.2$min.received.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Bittersweet"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Bittersweet"]) +
xlab('grooming rate subset 1') +
ylab('grooming rate subset 2') +
ggtitle('grooming received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep='')) #+
#stat_poly_eq(formula = y~x,
# aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
# parse = TRUE)
p1
p2
p3
Ok, looks pretty strongly correlated, now we’ll look at the agg behaviors
SGEI.intros <- intros
# aggression received
agg.function <- function(socialDF){
sge.groom.rate.df <- groom.df.1
aggs <- prox<-lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","TC","AT","CH")),]})
aggs.contact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TC","AT")),]})
aggs.nocontact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","CH")),]})
agg.received<-lapply(aggs,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.given<-lapply(aggs,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg.contact.received<-lapply(aggs.contact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.contact.given<-lapply(aggs.contact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg.nocontact.received<-lapply(aggs.nocontact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.nocontact.given<-lapply(aggs.nocontact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.received[[i]]$Group.1
tmp4 <- agg.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub <- SGEI.intros
intro_sub$agg.received <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub <- na.omit(intro_sub)
intro_sub$agg.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub$id.group))
tmp$agg.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.contact.received[[i]]$Group.1
tmp4 <- agg.contact.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub2 <- SGEI.intros
intro_sub2$agg.contact.received <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub2$id.group))
tmp$agg.contact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.nocontact.received[[i]]$Group.1
tmp4 <- agg.nocontact.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros
intro_sub3$agg.nocontact.received <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.received <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.received','num.obs')],intro_sub2$agg.contact.received,intro_sub3$agg.nocontact.received)
sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
tmp.final.behav.df <- sge.I.behavioral.df2[,c(4,6:7)]
# aggression given
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.given[[i]]$Group.1
tmp4 <- agg.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub <- SGEI.intros
intro_sub$agg.given <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub$agg.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub$id.group))
tmp$agg.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.contact.given[[i]]$Group.1
tmp4 <- agg.contact.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub2 <- SGEI.intros
intro_sub2$agg.contact.given <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub2$id.group))
tmp$agg.contact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros$group)){
tmp3 <- agg.nocontact.given[[i]]$Group.1
tmp4 <- agg.nocontact.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros
intro_sub3$agg.nocontact.given <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.given <- tmp$agg
# add zero rows back in
SGEI.intros$id.group <- paste(SGEI.intros$ID,SGEI.intros$group, sep='_')
tmp <- subset(SGEI.intros, !(SGEI.intros$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.given','num.obs')],intro_sub2$agg.contact.given,intro_sub3$agg.nocontact.given)
sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
final.behav.df <- cbind(sge.I.behavioral.df2,tmp.final.behav.df)
colnames(final.behav.df) <- c('ID',
'group',
'phase',
'agg.given',
'num.obs',
'agg.contact.given',
'agg.nocontact.given',
'agg.received',
'agg.contact.received',
'agg.nocontact.received')
final.behav.df$agg.given.rate <- (final.behav.df$agg.given / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.given.rate <- (final.behav.df$agg.contact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.given.rate <- (final.behav.df$agg.nocontact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.received.rate <- (final.behav.df$agg.received / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.received.rate <- (final.behav.df$agg.contact.received / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.received.rate <- (final.behav.df$agg.nocontact.received / final.behav.df$num.obs) * 2
return(final.behav.df)
}
agg.df.1 <- agg.function(social_1)
agg.df.2 <- agg.function(social_2)
tmp <- cor.test(agg.df.1$agg.given.rate,agg.df.2$agg.given.rate)
p1 <- ggplot(agg.df.1, aes(x=agg.given.rate,y=agg.df.2$agg.given.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg given rate subset 1') +
ylab('agg given rate subset 2') +
ggtitle('agg given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
tmp <- cor.test(agg.df.1$agg.contact.given.rate,agg.df.2$agg.contact.given.rate)
p2 <- ggplot(agg.df.1, aes(x=agg.contact.given.rate,y=agg.df.2$agg.contact.given.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg contact given rate subset 1') +
ylab('agg contact given rate subset 2') +
ggtitle('agg contact given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
tmp <- cor.test(agg.df.1$agg.nocontact.given.rate,agg.df.2$agg.nocontact.given.rate)
p3 <- ggplot(agg.df.1, aes(x=agg.nocontact.given.rate,y=agg.df.2$agg.nocontact.given.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg noncontact given rate subset 1') +
ylab('agg noncontact given rate subset 2') +
ggtitle('agg noncontact given rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
tmp <- cor.test(agg.df.1$agg.received.rate,agg.df.2$agg.received.rate)
p4 <- ggplot(agg.df.1, aes(x=agg.received.rate,y=agg.df.2$agg.received.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg received rate subset 1') +
ylab('agg received rate subset 2') +
ggtitle('agg received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
tmp <- cor.test(agg.df.1$agg.contact.received.rate,agg.df.2$agg.contact.received.rate)
p5 <- ggplot(agg.df.1, aes(x=agg.contact.received.rate,y=agg.df.2$agg.contact.received.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg contact received rate subset 1') +
ylab('agg contact received rate subset 2') +
ggtitle('agg contact received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
tmp <- cor.test(agg.df.1$agg.nocontact.received.rate,agg.df.2$agg.nocontact.received.rate)
p6 <- ggplot(agg.df.1, aes(x=agg.nocontact.received.rate,y=agg.df.2$agg.nocontact.received.rate)) +
geom_point(alpha = 0.6, color=brocolors("crayons")["Indigo"]) +
#geom_density_2d(color=brocolors("crayons")["Bittersweet"], size=1) +
geom_abline(slope = 1, intercept = 0, lty = 'dashed') +
stat_smooth(method='lm', color=brocolors("crayons")["Indigo"]) +
xlab('agg noncontact received rate subset 1') +
ylab('agg noncontact received rate subset 2') +
ggtitle('agg noncontact received rate',subtitle = paste('R = ',round(tmp$estimate,2),'; 95% CI: ',round(tmp$conf.int,2)[1],'-',round(tmp$conf.int,2)[2],'; p = ',round(tmp$p.value,20),sep=''))
p1
p2
p3
p4
p5
p6
Everything looks pretty well correlated. Perhaps as expected the noncontact agg measures are better correlated than contact agg measures.
Now let’s look at how the grooming and agg behaviors are distributed across a sliding window of obs. Here I’ve calculated grooming and agg data for 29 windows representing ~6 months of data collection and ~ 60 hours of data per group. The window size is 3 obs.
agg.function <- function(socialDF){
sge.groom.rate.df <- groom.df.1
SGEI.intros_tmp <- intros_tmp
aggs <- prox<-lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","TC","AT","CH")),]})
aggs.contact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TC","AT")),]})
#aggs.contact[[2]] <- aggs.contact[[1]]
#aggs.contact[[5]] <- aggs.contact[[1]]
#aggs.contact[[16]] <- aggs.contact[[1]]
for (i in 1:18){
if (nrow(aggs.contact[[i]]) == 0){
aggs.contact[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
}
}
for (i in 1:18){
if (nrow(aggs[[i]]) == 0){
aggs[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
}
}
aggs.nocontact <- lapply(socialDF,function(x){x[which(x[,5]%in%c("TN","CH")),]})
for (i in 1:18){
if (nrow(aggs.nocontact[[i]]) == 0){
aggs.nocontact[[i]][1,1:6] <- c('2013-10-10',1,1,1,1,1)
}
}
agg.received<-lapply(aggs,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.given<-lapply(aggs,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg.contact.received<-lapply(aggs.contact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.contact.given<-lapply(aggs.contact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg.nocontact.received<-lapply(aggs.nocontact,function(x){aggregate(x[,4],by=list(x[,4]),length)})
agg.nocontact.given<-lapply(aggs.nocontact,function(x){aggregate(x[,3],by=list(x[,3]),length)})
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.received[[i]]$Group.1
tmp4 <- agg.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub <- SGEI.intros_tmp
intro_sub$agg.received <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub <- na.omit(intro_sub)
intro_sub$agg.received <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub$id.group))
tmp$agg.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.contact.received[[i]]$Group.1
tmp4 <- agg.contact.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub2 <- SGEI.intros_tmp
intro_sub2$agg.contact.received <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.received <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub2$id.group))
tmp$agg.contact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.nocontact.received[[i]]$Group.1
tmp4 <- agg.nocontact.received[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros_tmp
intro_sub3$agg.nocontact.received <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.received <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.received <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.received','num.obs')],intro_sub2$agg.contact.received,intro_sub3$agg.nocontact.received)
sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
tmp.final.behav.df <- sge.I.behavioral.df2[,c(4,6:7)]
# aggression given
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.given[[i]]$Group.1
tmp4 <- agg.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub <- SGEI.intros_tmp
intro_sub$agg.given <- 0
intro_sub$id.group <- paste(intro_sub$ID,intro_sub$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub$id.group)
intro_sub <- subset(intro_sub, intro_sub$id.group %in% tmp$id.group)
intro_sub <- intro_sub[match(tmp$id.group, intro_sub$id.group),]
intro_sub$agg.given <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub$id.group))
tmp$agg.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub <- rbind.data.frame(intro_sub,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.contact.given[[i]]$Group.1
tmp4 <- agg.contact.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
intro_sub2 <- SGEI.intros_tmp
intro_sub2$agg.contact.given <- 0
intro_sub2$id.group <- paste(intro_sub2$ID,intro_sub2$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub2$id.group)
intro_sub2 <- subset(intro_sub2, intro_sub2$id.group %in% tmp$id.group)
intro_sub2 <- intro_sub2[match(tmp$id.group, intro_sub2$id.group),]
intro_sub2$agg.contact.given <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub2$id.group))
tmp$agg.contact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub2 <- rbind.data.frame(intro_sub2,tmp)
####
agg <- vector()
name <- vector()
group <- vector()
for (i in unique(SGEI.intros_tmp$group)){
tmp3 <- agg.nocontact.given[[i]]$Group.1
tmp4 <- agg.nocontact.given[[i]]$x
tmp5 <- rep(i,length(tmp3))
agg <- c(agg,tmp4)
name <- c(name,tmp3)
group <- c(group,tmp5)
}
tmp <- as.data.frame(cbind(agg,name,group))
tmp <- tmp[!(tmp$name=="ZZ"),]
tmp$agg <- as.numeric(as.character(tmp$agg))
tmp$name <- toupper(tmp$name)
tmp$id.group <- paste(tmp$name,tmp$group, sep='_')
####
intro_sub3 <- SGEI.intros_tmp
intro_sub3$agg.nocontact.given <- 0
intro_sub3$id.group <- paste(intro_sub3$ID,intro_sub3$group, sep='_')
tmp <- subset(tmp, tmp$id.group %in% intro_sub3$id.group)
intro_sub3 <- subset(intro_sub3, intro_sub3$id.group %in% tmp$id.group)
intro_sub3 <- intro_sub3[match(tmp$id.group, intro_sub3$id.group),]
intro_sub3$agg.nocontact.given <- tmp$agg
# add zero rows back in
SGEI.intros_tmp$id.group <- paste(SGEI.intros_tmp$ID,SGEI.intros_tmp$group, sep='_')
tmp <- subset(SGEI.intros_tmp, !(SGEI.intros_tmp$id.group %in% intro_sub3$id.group))
tmp$agg.nocontact.given <- 0
#tmp <- tmp[,c(1:8,10,9)]
intro_sub3 <- rbind.data.frame(intro_sub3,tmp)
intro_sub <- intro_sub[order(intro_sub$id.group),]
intro_sub2 <- intro_sub2[order(intro_sub2$id.group),]
intro_sub3 <- intro_sub3[order(intro_sub3$id.group),]
sge.I.behavioral.df2 <- cbind(intro_sub[,c('ID','group','phase','agg.given','num.obs')],intro_sub2$agg.contact.given,intro_sub3$agg.nocontact.given)
sge.I.behavioral.df2 <- sge.I.behavioral.df2[order(sge.I.behavioral.df2$phase,sge.I.behavioral.df2$group,sge.I.behavioral.df2$ID),]
sge.groom.rate.df <- sge.groom.rate.df[order(sge.groom.rate.df$phase,sge.groom.rate.df$group,sge.groom.rate.df$ID),]
final.behav.df <- cbind(sge.I.behavioral.df2,tmp.final.behav.df)
colnames(final.behav.df) <- c('ID',
'group',
'phase',
'agg.given',
'num.obs',
'agg.contact.given',
'agg.nocontact.given',
'agg.received',
'agg.contact.received',
'agg.nocontact.received')
final.behav.df$agg.given.rate <- (final.behav.df$agg.given / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.given.rate <- (final.behav.df$agg.contact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.given.rate <- (final.behav.df$agg.nocontact.given / final.behav.df$num.obs) * 2
final.behav.df$agg.received.rate <- (final.behav.df$agg.received / final.behav.df$num.obs) * 2
final.behav.df$agg.contact.received.rate <- (final.behav.df$agg.contact.received / final.behav.df$num.obs) * 2
final.behav.df$agg.nocontact.received.rate <- (final.behav.df$agg.nocontact.received / final.behav.df$num.obs) * 2
return(final.behav.df)
}
#sliding window
windows.matrix <- matrix(nrow = 27, ncol = 3)
for (i in 1:27){
windows.matrix[i,] <- seq(i,i+2,1)
}
windows.matrix <- cbind(windows.matrix, paste(windows.matrix[,1],':',windows.matrix[,3],sep=''))
for (i in 1:27){
tmp <- lapply(names(data_row),function(x){
tmpRange <- windows.matrix[i,4]
tmp2 <- subset(data_row[[x]], obs.no %in% seq(i,i+2,1))
return(tmp2)
})
names(tmp) <- names(data_row)
assign(value = tmp, x = paste('data_window_',i,sep=''))
}
groom.function2 <- function(groomDF){
groom.tot1<-data.frame(matrix(ncol=19,nrow=360))
colnames(groom.tot1)<-c("ID1","ID2","group","no.obs","bouts","ID1.groom","ID2.groom","min","ID1.groom.min","ID2.groom.min","ID1.rank","ID2.rank","proximity","contact.min","CSI","elo.ID1","elo.ID2","ID1.agg","ID2.agg")
groom.tot1[,"ID1"]<-rep(intros_tmp[,"ID"],4)
for (i in 1:180) {groom.tot1[i,"group"]<-intros_tmp[which(intros_tmp[1:45,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
for (i in 181:360) {groom.tot1[i,"group"]<-intros_tmp[46:90,][which(intros_tmp[46:90,"ID"]==groom.tot1[i,1]),"group"]}; rm(i)
groom.tot1<-groom.tot1[order(groom.tot1[,3],groom.tot1[,1]),];rownames(groom.tot1)<-NULL
for (i in unique(intros_tmp[,4])) {intros_tmp[1:45,][which(intros_tmp[1:45,"ID"]==i),"group"]->grp;
groom.tot1[grep("SGE1",groom.tot1$group),][which(groom.tot1[grep("SGE1",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros_tmp[1:45,][which(intros_tmp[1:45,"group"]==grp),"ID"],i)}
for (i in unique(intros_tmp[,4])) {intros_tmp[46:90,][which(intros_tmp[46:90,"ID"]==i),"group"]->grp;
groom.tot1[grep("SGE2",groom.tot1$group),][which(groom.tot1[grep("SGE2",groom.tot1$group),"ID1"]==i),2]<-setdiff(intros_tmp[46:90,][which(intros_tmp[46:90,"group"]==grp),"ID"],i)}
groom.tot<-groom.tot1; rm(groom.tot1)
groom.tot[,1]<-as.character(groom.tot[,1]); groom.tot[,2]<-as.character(groom.tot[,2]); groom.tot[,3]<-as.character(groom.tot[,3]); rm(grp); rm(i)
##caclulate grooming bouts
groom.tot$ID1.intro.no<-NA; groom.tot$ID2.intro.no<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
b1=sum(g[,5]=="GM"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="GM"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,5]=b1+b2;groom.tot[id,6]=b1;groom.tot[id,7]=b2
groom.tot[id,4]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"num.obs"]
groom.tot[id,"ID1.rank"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"finalElo"]
groom.tot[id,"ID2.rank"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"finalElo"]
groom.tot[id,"ID1.intro.no"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"intro.no"]
groom.tot[id,"ID2.intro.no"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"intro.no"]
}#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
groom.tot$ID1.intro.date<-NA; groom.tot$ID2.intro.date<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0;grp<-groom.tot[id,"group"]
groom.tot[id,"ID1.intro.date"]=intros_tmp[which(intros_tmp[,"ID"]==id1&intros_tmp$group==grp),"intro.date"]
groom.tot[id,"ID2.intro.date"]=intros_tmp[which(intros_tmp[,"ID"]==id2&intros_tmp$group==grp),"intro.date"]
}#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
## calculate who ended the grooming bout
groom.tot$ID1.end.bout<-NA; groom.tot$ID2.end.bout<-NA
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-groomDF[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
b1=sum(g[,5]=="G-"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="G-"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,"ID1.end.bout"]=b1
groom.tot[id,"ID2.end.bout"]=b2
}#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
groom.tot$ID1.approach<-NA
groom.tot$ID2.approach<-NA
##caclulate approaches
for (id in 1:nrow(groom.tot)){
id1<-groom.tot[id,"ID1"]; id2<-groom.tot[id,"ID2"]; g<-prox[[as.character(groom.tot[id,"group"])]]; b1=0;b2=0
b1=sum(g[,5]=="PX"&g[,3]==id1 & g[,4]==id2)
b2=sum(g[,5]=="PX"&g[,4]==id1 & g[,3]==id2)
groom.tot[id,"ID1.approach"]=b1;groom.tot[id,"ID2.approach"]=b2
}#; rm(g);rm(b1);rm(b2);rm(id1);rm(id2);rm(id)
## calculate grooming times
groom.tot$ID1.groom.min<-NA
groom.tot$ID2.groom.min<-NA
groom.tot$min<-NA
for (i in 1:nrow(groom.tot)){
id1<-groom.tot[i,"ID1"]; id2<-groom.tot[i,"ID2"]; t1=0;t2=0
g<-groomDF[[as.character(groom.tot[i,"group"])]]
if (id1%in%c(g[,3],g[,4])&id2%in%c(g[,3],g[,4])){
g<-g[which((g[,4]==id1 & g[,3]==id2 )|(g[,3]==id1 & g[,4]==id2)|g[,3]=="ZZ"),]
for (r in 1:nrow(g)){
if (g[r,5]=="GM" & g[r,3]==id1 & g[r,4]==id2) t1=t1+(g[r+1,2]-g[r,2])
if (g[r,5]=="GM" & g[r,4]==id1 & g[r,3]==id2) t2=t2+(g[r+1,2]-g[r,2])
}
groom.tot[i,"ID1.groom.min"]=t1
groom.tot[i,"ID2.groom.min"]=t2
groom.tot[i,"min"]=t1+t2}
}
#rm(t1);rm(t2);rm(i);rm(id1);rm(id2);rm(r);rm(g)
###########
phase <- groom.tot$group
phase <- substr(phase,(nchar(phase)+1)-4,nchar(phase))
groom.tot$phase <- phase
## generate a data frame with minutes of grooming (total), elo, # of obs, and group, per individual
## phase I
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
tmp <- subset(groom.tot, groom.tot$phase == 'SGE1')
tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
tmp.minutes <- cbind(tmp.minutes,tmp2)
tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
tmp.names <- cbind(tmp.names,i)
tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}
sge.I.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.I.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.I.groom.rate.df$elo <- as.numeric(as.character(sge.I.groom.rate.df$elo))
sge.I.groom.rate.df$min <- as.numeric(as.character(sge.I.groom.rate.df$min))
sge.I.groom.rate.df$min.given <- as.numeric(as.character(sge.I.groom.rate.df$min.given))
sge.I.groom.rate.df$min.received <- as.numeric(as.character(sge.I.groom.rate.df$min.received))
sge.I.groom.rate.df$no.obs <- as.numeric(as.character(sge.I.groom.rate.df$no.obs))
sge.I.groom.rate.df$phase <- 1
## now for phase II
tmp.minutes <- vector()
tmp.minutes.given <- vector()
tmp.minutes.received <- vector()
tmp.names <- vector()
tmp.ranks <- vector()
tmp.obs <- vector()
tmp.group <- vector()
for (i in unique(groom.tot$ID1)){
tmp <- subset(groom.tot, groom.tot$phase == 'SGE2')
tmp2 <- sum(subset(tmp, tmp$ID1 == i)$min)
tmp3 <- sum(subset(tmp, tmp$ID1 == i)$ID1.groom.min)
tmp4 <- sum(subset(tmp, tmp$ID1 == i)$ID2.groom.min)
tmp.minutes <- cbind(tmp.minutes,tmp2)
tmp.minutes.given <- cbind(tmp.minutes.given,tmp3)
tmp.minutes.received <- cbind(tmp.minutes.received,tmp4)
tmp.names <- cbind(tmp.names,i)
tmp.ranks <- cbind(tmp.ranks,subset(tmp, tmp$ID1 == i)$ID1.rank[1])
tmp.obs <- cbind(tmp.obs,subset(tmp, tmp$ID1 == i)$no.obs[1])
tmp.group <- cbind(tmp.group,subset(tmp, tmp$ID1 == i)$group[1])
}
sge.II.groom.rate.df <- as.data.frame(t(rbind(tmp.names,tmp.ranks,tmp.minutes,tmp.minutes.given,tmp.minutes.received,tmp.obs,tmp.group)))
names(sge.II.groom.rate.df) <- c('ID','elo','min','min.given','min.received','no.obs','group')
sge.II.groom.rate.df$elo <- as.numeric(as.character(sge.II.groom.rate.df$elo))
sge.II.groom.rate.df$min <- as.numeric(as.character(sge.II.groom.rate.df$min))
sge.II.groom.rate.df$min.given <- as.numeric(as.character(sge.II.groom.rate.df$min.given))
sge.II.groom.rate.df$min.received <- as.numeric(as.character(sge.II.groom.rate.df$min.received))
sge.II.groom.rate.df$no.obs <- as.numeric(as.character(sge.II.groom.rate.df$no.obs))
sge.II.groom.rate.df$phase <- 2
sge.groom.rate.df <- rbind(sge.I.groom.rate.df,sge.II.groom.rate.df)
return(sge.groom.rate.df)
}
for (i in ls(pattern='data_window')){
tmp <- get(i)
intros_tmp <- intros
for (x in names(data)) intros_tmp[which(intros_tmp[,"group"]==x),"num.obs"]=length(levels(as.factor((tmp[[x]][,"obs.no"]))))
social_tmp <- lapply(tmp,convert.soc)
groom_tmp <- lapply(social_tmp,function(x){x[which(x[,5]%in%c("GM","G-","ZZ")),]})
names(groom_tmp) <- names(tmp)
groom_tmp2 <- groom.function2(groom_tmp)
agg_tmp <- agg.function(social_tmp)
assign(value = groom_tmp2, x = paste('groom_df_window_',strsplit(i, split='_')[[1]][3],sep=''))
assign(value = agg_tmp, x = paste('agg_df_window_',strsplit(i, split='_')[[1]][3],sep=''))
}
test_minute_matrix <- matrix(nrow = 90, ncol = 27)
#Pattern1 <- grep("df_window",names(.GlobalEnv),value=TRUE)
test_list <- list(groom_df_window_1,groom_df_window_2,groom_df_window_3,groom_df_window_4,groom_df_window_5,groom_df_window_6,groom_df_window_7,groom_df_window_8,groom_df_window_9,groom_df_window_10,groom_df_window_11,groom_df_window_12,groom_df_window_13,groom_df_window_14,groom_df_window_15,groom_df_window_16,groom_df_window_17,groom_df_window_18,groom_df_window_19,groom_df_window_20,groom_df_window_21,groom_df_window_22,groom_df_window_23,groom_df_window_24,groom_df_window_25,groom_df_window_26,groom_df_window_27)
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,3]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'minutes grooming')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,4]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'minutes grooming given')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,5]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'minutes grooming received')
# agg plots
test_minute_matrix <- matrix(nrow = 90, ncol = 27)
test_list <- list(agg_df_window_1,agg_df_window_2,agg_df_window_3,agg_df_window_4,agg_df_window_5,agg_df_window_6,agg_df_window_7,agg_df_window_8,agg_df_window_9,agg_df_window_10,agg_df_window_11,agg_df_window_12,agg_df_window_13,agg_df_window_14,agg_df_window_15,agg_df_window_16,agg_df_window_17,agg_df_window_18,agg_df_window_19,agg_df_window_20,agg_df_window_21,agg_df_window_22,agg_df_window_23,agg_df_window_24,agg_df_window_25,agg_df_window_26,agg_df_window_27)
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,4]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'aggs given')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,8]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'aggs received')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,6]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'contact aggs given')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,7]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'noncontact aggs given')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,9]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'contact aggs received')
for (i in 1:27){
test_minute_matrix[,i] <- test_list[[i]][,10]
}
test_minute_matrix_melt <- reshape2::melt(test_minute_matrix)
ggplot(test_minute_matrix_melt, aes(x=as.factor(Var2),y=value)) +
geom_boxplot() +
xlab('sliding window') +
ggtitle(label='',subtitle = 'noncontact aggs received')
#### Produce heatmap of FDR sig. genes from all behavioral models on HARDAC
results_matrix_tmp <- matrix(nrow = 13, ncol = 14)
rownames(results_matrix_tmp) <- c('elo',
'aggDiff',
'Contact agg given',
'Contact agg rec',
'Groom given',
'Groom rec',
'Noncontact agg given',
'Noncontact agg rec',
'Overall agg given',
'Overall agg rec',
'Overall groom',
'PC1.agg',
'PC1.groom')
colnames(results_matrix_tmp) <- c("CD14",
"CD16",
"CD20",
"CD4",
"CD8",
"dexCA",
"dexCA_delta",
"dexGE",
"dexGE_delta",
"gard",
"gard_delta",
"lps",
"lps_delta",
"null")
elo_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_FDR.rds')
elo_CD4 <- elo_flow[[1]]
elo_CD8 <- elo_flow[[2]]
elo_CD14 <- elo_flow[[3]]
elo_CD16 <- elo_flow[[4]]
elo_CD20 <- elo_flow[[5]]
rm(elo_flow)
elo_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_gard_modelRes_fdr')
elo_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/Elo_GARDdelta_modelRes.txt')
elo_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_Elo.txt')
elo_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/Elo_LPSdelta_modelRes.txt')
elo_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_ca_modelRes_fdr')
elo_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_CA_delta_modelRes_fdr')
elo_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_ge_modelRes_fdr')
elo_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/elo_GE_delta_modelRes_fdr')
PC1.agg_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_FDR.rds')
PC1.agg_CD4 <- PC1.agg_flow[[1]]
PC1.agg_CD8 <- PC1.agg_flow[[2]]
PC1.agg_CD14 <- PC1.agg_flow[[3]]
PC1.agg_CD16 <- PC1.agg_flow[[4]]
PC1.agg_CD20 <- PC1.agg_flow[[5]]
rm(PC1.agg_flow)
PC1.agg_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_gard_modelRes_fdr')
PC1.agg_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_GARDdelta_modelRes.txt')
PC1.agg_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_PC1.agg.txt')
PC1.agg_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_LPSdelta_modelRes.txt')
PC1.agg_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_ca_modelRes_fdr')
PC1.agg_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_CA_delta_modelRes_fdr')
PC1.agg_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_ge_modelRes_fdr')
PC1.agg_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.agg_GE_delta_modelRes_fdr')
PC1.groom_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_FDR.rds')
PC1.groom_CD4 <- PC1.groom_flow[[1]]
PC1.groom_CD8 <- PC1.groom_flow[[2]]
PC1.groom_CD14 <- PC1.groom_flow[[3]]
PC1.groom_CD16 <- PC1.groom_flow[[4]]
PC1.groom_CD20 <- PC1.groom_flow[[5]]
rm(PC1.groom_flow)
PC1.groom_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_gard_modelRes_fdr')
PC1.groom_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_GARDdelta_modelRes.txt')
PC1.groom_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_PC1.groom.txt')
PC1.groom_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_LPSdelta_modelRes.txt')
PC1.groom_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_ca_modelRes_fdr')
PC1.groom_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_CA_delta_modelRes_fdr')
PC1.groom_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_ge_modelRes_fdr')
PC1.groom_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/PC1.groom_GE_delta_modelRes_fdr')
gc.agg.Diff_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_FDR.rds')
gc.agg.Diff_CD4 <- gc.agg.Diff_flow[[1]]
gc.agg.Diff_CD8 <- gc.agg.Diff_flow[[2]]
gc.agg.Diff_CD14 <- gc.agg.Diff_flow[[3]]
gc.agg.Diff_CD16 <- gc.agg.Diff_flow[[4]]
gc.agg.Diff_CD20 <- gc.agg.Diff_flow[[5]]
rm(gc.agg.Diff_flow)
gc.agg.Diff_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_gard_modelRes_fdr')
gc.agg.Diff_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_GARDdelta_modelRes.txt')
gc.agg.Diff_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.agg.Diff.txt')
gc.agg.Diff_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_LPSdelta_modelRes.txt')
gc.agg.Diff_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_ca_modelRes_fdr')
gc.agg.Diff_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_CA_delta_modelRes_fdr')
gc.agg.Diff_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_ge_modelRes_fdr')
gc.agg.Diff_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.agg.Diff_GE_delta_modelRes_fdr')
gc.contact.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_FDR.rds')
gc.contact.agg.given_CD4 <- gc.contact.agg.given_flow[[1]]
gc.contact.agg.given_CD8 <- gc.contact.agg.given_flow[[2]]
gc.contact.agg.given_CD14 <- gc.contact.agg.given_flow[[3]]
gc.contact.agg.given_CD16 <- gc.contact.agg.given_flow[[4]]
gc.contact.agg.given_CD20 <- gc.contact.agg.given_flow[[5]]
rm(gc.contact.agg.given_flow)
gc.contact.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_gard_modelRes_fdr')
gc.contact.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_GARDdelta_modelRes.txt')
gc.contact.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.contact.agg.given.txt')
gc.contact.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_LPSdelta_modelRes.txt')
gc.contact.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_ca_modelRes_fdr')
gc.contact.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_CA_delta_modelRes_fdr')
gc.contact.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_ge_modelRes_fdr')
gc.contact.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.given_GE_delta_modelRes_fdr')
gc.contact.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_FDR.rds')
gc.contact.agg.rec_CD4 <- gc.contact.agg.rec_flow[[1]]
gc.contact.agg.rec_CD8 <- gc.contact.agg.rec_flow[[2]]
gc.contact.agg.rec_CD14 <- gc.contact.agg.rec_flow[[3]]
gc.contact.agg.rec_CD16 <- gc.contact.agg.rec_flow[[4]]
gc.contact.agg.rec_CD20 <- gc.contact.agg.rec_flow[[5]]
rm(gc.contact.agg.rec_flow)
gc.contact.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_gard_modelRes_fdr')
gc.contact.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_GARDdelta_modelRes.txt')
gc.contact.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.contact.agg.rec.txt')
gc.contact.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_LPSdelta_modelRes.txt')
gc.contact.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_ca_modelRes_fdr')
gc.contact.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_CA_delta_modelRes_fdr')
gc.contact.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_ge_modelRes_fdr')
gc.contact.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.contact.agg.rec_GE_delta_modelRes_fdr')
gc.groom.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_FDR.rds')
gc.groom.given_CD4 <- gc.groom.given_flow[[1]]
gc.groom.given_CD8 <- gc.groom.given_flow[[2]]
gc.groom.given_CD14 <- gc.groom.given_flow[[3]]
gc.groom.given_CD16 <- gc.groom.given_flow[[4]]
gc.groom.given_CD20 <- gc.groom.given_flow[[5]]
rm(gc.groom.given_flow)
gc.groom.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_gard_modelRes_fdr')
gc.groom.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_GARDdelta_modelRes.txt')
gc.groom.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.groom.given.txt')
gc.groom.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_LPSdelta_modelRes.txt')
gc.groom.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_ca_modelRes_fdr')
gc.groom.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_CA_delta_modelRes_fdr')
gc.groom.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_ge_modelRes_fdr')
gc.groom.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.given_GE_delta_modelRes_fdr')
gc.groom.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_FDR.rds')
gc.groom.rec_CD4 <- gc.groom.rec_flow[[1]]
gc.groom.rec_CD8 <- gc.groom.rec_flow[[2]]
gc.groom.rec_CD14 <- gc.groom.rec_flow[[3]]
gc.groom.rec_CD16 <- gc.groom.rec_flow[[4]]
gc.groom.rec_CD20 <- gc.groom.rec_flow[[5]]
rm(gc.groom.rec_flow)
gc.groom.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_gard_modelRes_fdr')
gc.groom.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_GARDdelta_modelRes.txt')
gc.groom.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.groom.rec.txt')
gc.groom.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_LPSdelta_modelRes.txt')
gc.groom.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_ca_modelRes_fdr')
gc.groom.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_CA_delta_modelRes_fdr')
gc.groom.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_ge_modelRes_fdr')
gc.groom.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.groom.rec_GE_delta_modelRes_fdr')
gc.noncontact.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_FDR.rds')
gc.noncontact.agg.given_CD4 <- gc.noncontact.agg.given_flow[[1]]
gc.noncontact.agg.given_CD8 <- gc.noncontact.agg.given_flow[[2]]
gc.noncontact.agg.given_CD14 <- gc.noncontact.agg.given_flow[[3]]
gc.noncontact.agg.given_CD16 <- gc.noncontact.agg.given_flow[[4]]
gc.noncontact.agg.given_CD20 <- gc.noncontact.agg.given_flow[[5]]
rm(gc.noncontact.agg.given_flow)
gc.noncontact.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_gard_modelRes_fdr')
gc.noncontact.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_GARDdelta_modelRes.txt')
gc.noncontact.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.noncontact.agg.given.txt')
gc.noncontact.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_LPSdelta_modelRes.txt')
gc.noncontact.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_ca_modelRes_fdr')
gc.noncontact.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_CA_delta_modelRes_fdr')
gc.noncontact.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_ge_modelRes_fdr')
gc.noncontact.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.given_GE_delta_modelRes_fdr')
gc.noncontact.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_FDR.rds')
gc.noncontact.agg.rec_CD4 <- gc.noncontact.agg.rec_flow[[1]]
gc.noncontact.agg.rec_CD8 <- gc.noncontact.agg.rec_flow[[2]]
gc.noncontact.agg.rec_CD14 <- gc.noncontact.agg.rec_flow[[3]]
gc.noncontact.agg.rec_CD16 <- gc.noncontact.agg.rec_flow[[4]]
gc.noncontact.agg.rec_CD20 <- gc.noncontact.agg.rec_flow[[5]]
rm(gc.noncontact.agg.rec_flow)
gc.noncontact.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_gard_modelRes_fdr')
gc.noncontact.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_GARDdelta_modelRes.txt')
gc.noncontact.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.noncontact.agg.rec.txt')
gc.noncontact.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_LPSdelta_modelRes.txt')
gc.noncontact.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_ca_modelRes_fdr')
gc.noncontact.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_CA_delta_modelRes_fdr')
gc.noncontact.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_ge_modelRes_fdr')
gc.noncontact.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.noncontact.agg.rec_GE_delta_modelRes_fdr')
gc.overall.agg.rec_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_FDR.rds')
gc.overall.agg.rec_CD4 <- gc.overall.agg.rec_flow[[1]]
gc.overall.agg.rec_CD8 <- gc.overall.agg.rec_flow[[2]]
gc.overall.agg.rec_CD14 <- gc.overall.agg.rec_flow[[3]]
gc.overall.agg.rec_CD16 <- gc.overall.agg.rec_flow[[4]]
gc.overall.agg.rec_CD20 <- gc.overall.agg.rec_flow[[5]]
rm(gc.overall.agg.rec_flow)
gc.overall.agg.rec_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_gard_modelRes_fdr')
gc.overall.agg.rec_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_GARDdelta_modelRes.txt')
gc.overall.agg.rec_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.agg.rec.txt')
gc.overall.agg.rec_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_LPSdelta_modelRes.txt')
gc.overall.agg.rec_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_ca_modelRes_fdr')
gc.overall.agg.rec_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_CA_delta_modelRes_fdr')
gc.overall.agg.rec_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_ge_modelRes_fdr')
gc.overall.agg.rec_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.rec_GE_delta_modelRes_fdr')
gc.overall.agg.given_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_FDR.rds')
gc.overall.agg.given_CD4 <- gc.overall.agg.given_flow[[1]]
gc.overall.agg.given_CD8 <- gc.overall.agg.given_flow[[2]]
gc.overall.agg.given_CD14 <- gc.overall.agg.given_flow[[3]]
gc.overall.agg.given_CD16 <- gc.overall.agg.given_flow[[4]]
gc.overall.agg.given_CD20 <- gc.overall.agg.given_flow[[5]]
rm(gc.overall.agg.given_flow)
gc.overall.agg.given_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_gard_modelRes_fdr')
gc.overall.agg.given_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_GARDdelta_modelRes.txt')
gc.overall.agg.given_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.agg.given.txt')
gc.overall.agg.given_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_LPSdelta_modelRes.txt')
gc.overall.agg.given_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_ca_modelRes_fdr')
gc.overall.agg.given_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_CA_delta_modelRes_fdr')
gc.overall.agg.given_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_ge_modelRes_fdr')
gc.overall.agg.given_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.agg.given_GE_delta_modelRes_fdr')
gc.overall.groom_flow <- readRDS('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_FDR.rds')
gc.overall.groom_CD4 <- gc.overall.groom_flow[[1]]
gc.overall.groom_CD8 <- gc.overall.groom_flow[[2]]
gc.overall.groom_CD14 <- gc.overall.groom_flow[[3]]
gc.overall.groom_CD16 <- gc.overall.groom_flow[[4]]
gc.overall.groom_CD20 <- gc.overall.groom_flow[[5]]
rm(gc.overall.groom_flow)
gc.overall.groom_gard <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_gard_modelRes_fdr')
gc.overall.groom_gard_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_GARDdelta_modelRes.txt')
gc.overall.groom_lps <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/lps_gc.overall.groom.txt')
gc.overall.groom_lps_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_LPSdelta_modelRes.txt')
gc.overall.groom_dexCA <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_ca_modelRes_fdr')
gc.overall.groom_dexCA_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_CA_delta_modelRes_fdr')
gc.overall.groom_dexGE <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_ge_modelRes_fdr')
gc.overall.groom_dexGE_delta <- read.table('/Users/noah/Desktop/SGE_behavioral_FDRmodelRes_tables/gc.overall.groom_GE_delta_modelRes_fdr')
# template code for inputting # sig genes
#countFlag <- 1
#for (i in ls(pattern = 'elo_')){
# tmp <- get(i)
# results_matrix_tmp[1,countFlag] <- sum(tmp[,ncol(tmp)] < 0.1)
# countFlag <- countFlag + 1
#}
#results_matrix_tmp[1,14] <- sum(elo_lps$fdr_Elo_at_NC < 0.1)
behaviorNames <- c("elo_",
"gc.agg.Diff_",
"gc.contact.agg.given_",
"gc.contact.agg.rec_",
"gc.groom.given_",
"gc.groom.rec_",
"gc.noncontact.agg.given_",
"gc.noncontact.agg.rec_",
"gc.overall.agg.given_",
"gc.overall.agg.rec_",
"gc.overall.groom_",
"PC1.agg_",
"PC1.groom_")
countFlag_beh <- 1
countFlag_dataset <- 1
for (b in behaviorNames){
tmp_filelist <- ls(pattern=b)
countFlag_dataset <- 1
for (i in tmp_filelist){
tmp <- get(i)
results_matrix_tmp[countFlag_beh,countFlag_dataset] <- sum(tmp[,ncol(tmp)] < 0.1)
countFlag_dataset <- countFlag_dataset + 1
}
countFlag_beh <- countFlag_beh + 1
}
countFlag_beh <- 1
for (i in behaviorNames){
tmp <- paste0(i, 'lps', sep='')
tmp <- get(tmp)
results_matrix_tmp[countFlag_beh,14] <- sum(tmp[,7] < 0.1)
countFlag_beh <- countFlag_beh + 1
}
results_matrix <- results_matrix_tmp[2:13,]
for (i in 1:ncol(results_matrix)){
results_matrix[,i] <- log2((results_matrix[,i] +1) / (results_matrix_tmp[1,i] +1))
}
colnames(results_matrix) <- c("CD14 flow-sorted (0)",
"CD16 flow-sorted (1635)",
"CD20 flow-sorted (83)",
"CD4 flow-sorted (326)",
"CD8 flow-sorted (23)",
"Dex CA (1)",
"Dex CA_delta (1)",
"Dex GE (1)",
"Dex GE_delta (1)",
"Gard (3517)",
"Gard_delta (0)",
"LPS (6234)",
"LPS_delta (4957)",
"NULL (4487)")
pheatmap(as.matrix(t(results_matrix)), color=colorRampPalette(c("lightblue4", "white", 'mediumpurple4'))(100), border_color = 'white', scale = 'row', legend = T, treeheight_row = 0, treeheight_col = 0, main = 'log2(FC) # of sig. genes over Elo', cluster_rows = 0, cluster_cols = 1)